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Abstract 

Three-dimensional simulations of the off-axis 
collisions of two drops are presented. The full 
Navier-Stokes equations are solved by a Front- 
Tracking/Finite-Difference method that allows a 
fully deformable fluid interface and the inclusion of 
surface tension. The drops are accelerated towards 
each other by a body force that is turned off before 
the drops collide. Depending on whether the 
interface between the drops is ruptured or not, the 
drops either bounce or coalesce. For drops that 
coalesce, the impact parameter, which measures 
how far the drops are off the symmetry line, 
determines the eventual outcome of the collision. 
For low impact parameters, the drops coalesce 
permanently, but for higher impact parameters, a 
grazing collision, where the drops coalesce and then 
stretch apart again is observed. The results are in 
agreement with experimental observations. 


Introduction 

The dynamic of fluid drops is of considerable 
importance in a number of engineering applications 
and natural processes. The combustion of fuel 
sprays, spray painting, various coating processes, 
as well as rain, are only a few of the more common 
examples. While it is usually the collective 
behavior of many drops that is of interest, often it 
is the motion of individual drops that determines 
the large scale properties of the system. Thus, for 
example, the total surface area of sprays depends on 
the size of the individual drops as well as their 
number density. Computational models for 
engineering predictions of spray combustion 
generally do not resolve the motion of individual 
drops and must rely on "subgrid" models where the 
average effects of the unresolved scales are 
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incorporated into the equations used to predict the 
large scale behavior. Many spray models (see 
Heywood 1 , for a discussion and references) use 
point particles to represent the drops. The drop 
motion is related to the fluid flow by empirical 
laws for drag, heat transfer and combustion. Often 
it is possible to focus on the dynamic of a single 
drop and how it interacts with the surrounding 
flow. When the number of drops per unit volume is 
high, however, it is necessary to account for the 
interactions between the drops and their collective 
effect on the flow. To account for drop collisions, 
models must contain “collision rules” that 
determine whether the drops coalesce or not. These 
rules are usually based on experimental 
investigations of binary collisions of drops, but the 
small spatial and temporal scales make detailed 
experimental measurements difficult and usually the 
record consist of little more than photographs or a 
video tape. Since the collision process generally 
involves large drop deformation and rupture of the 
interface separating the drops, it has not been 
amenable to detailed theoretical analysis. Previous 
studies are therefore mostly experimental, but 
sometimes supplemented by greatly simplified 
theoretical argument. 

Two recent experimental investigations of drop 
collisions can be found in Azhgriz and Poo 2 , and 
Jiang, Umemura and Law 3 who show several 
photographs of the various collision modes for 
both water and hydrocarbon drops. These, and other 
experimental investigations have provided 
considerable information and, in particular, it is 
now understood that the outcome of a collision can 
be classified into about five main categories. For 
head-on collisions we have four main categories: 
bouncing collision, where the drops collide and 
separate, retaining their identity; coalescence 
collision, where two drops become one; separation 
collision, where the drops temporarily become one 
but then break up again; and shattering collision, 
where the impact is so strong that the drops break 
up into several smaller drops. These categories 
survive for off-axis collisions, but a fifth one. 
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grazing or stretching collision, appears. Here, the 
drops coalesce upon contact, but are sufficiently far 
apart so that they continue along the original path 
and separate again. The form of the collision 
depends on the size of the drops, their relative 
velocities, their off-axis position and the physical 
properties of the fluids involved. For a given fluid, 
some of these collision regimes are not observer. 
Water drops, for example, do not show bouncing. 
(Jiang et a/. 3 , state that they also did not find 
reflective collision for water drops. This is 
apparently due to a limited parameter range studied 
by them as the experiments by Azhgriz and Poo 2 , 
show.) Other investigations of drop collisions may 
be found in Bradley and Stow 4 , and Podvysotsky 
and Shraiber 5 , for example. The major goal of these 
investigations has been to clarify the boundaries 
between the major collision categories and explain 
how they depend on the parameters of the problem. 
Simple models used to rationalize experimental 
findings have been presented by Park and Blair 6 , 
Ryley and Bennett-Cowell 7 , Brazier-Smith et al. 8 , 
Azhgriz and Poo 2 , and Jiang, Umemura and Law 3 . 

In principle, numerical solutions of the Navier- 
Strokes equations, where all scales of motion are 
fully resolved, can provide the missing 
information, but various numerical difficulties 
associated with moving boundaries between two 
fluids have made detailed simulations difficult in 
the past. Nevertheless, several authors have 
computed the axisymmetric head-on collision of 
drops with a wall. The earliest work is Foote 9 who 
followed the evolution of a rebounding 
axisymmetric drops at low Weber number using the 
MAC method. More recent computations work can 
be found in Fukai et al 10 who use a moving finite 
element method. We have recently conducted a 
numerical study of the head-on collision of two 
axisymmetric drops, see Nobari, Jan and 
Tryggvason 11 , where we examined the boundary 
between coalescing and reflecting collision for equal 
size drops. Here, we present numerical simulations 
of three-dimensional, off-axis collisions, where the 
full Navier Stokes equations are solved to give a 
detailed picture of the flow during collision. 

Formulation and Numerical Method 

The numerical technique used for the simulations 
presented in this paper is the Front Tracking/Finite 
Difference method of Unverdi and Tryggvason 12 * 13 . 
Since the procedure has bee described in detail 
before, we only outline it briefly here. 

The physical problem and the computational 
domain is sketched in Figure 1. The domain is a 


rectangular box and the drops are initially placed 
near each end of the domain. A force that is turned 
off before the drops collide, is applied to drive them 
together initially. Generally, the density and 
viscosity of the ambient fluid are much smaller 
than of the drop fluid and thus have only a small 
effect on the results. While it is therefore often 
sufficient to solve only for the fluid motion inside 
the drop, here we solve for the motion everywhere, 
both inside and outside the drops. The Navier- 
Stokes equations are valid for both fluids, and a 
single set of equations can be written for the whole 
domain as long as the jump in viscosity and 
density is correctly accounted for and surface 
tension is included: 

V -puu = -Vp + f x 

+ V ■ p| ViT + Vw T j + F a - Xf ). 

Here, u is the velocity, p is the pressure, and p 
and p are the discontinuous density and viscosity 
fields, respectively. F a is the surface tension force 
and f x is a body force used to give the drops their 
initial velocity. Notice that the surface tension 
force has been added as a delta function, only 
affecting the equations where the interface is. The 
detailed form of F 0 will be discussed below. The 
above equations are supplemented by the 
incompressibility conditions 
V-w =0 

which, when combined with the momentum 
equations leads to a non-separable elliptic equation 
for the pressure. We also have equations of state for 
the density and viscosity: 

&+u-Vp = 0 

Ot 

iE+u -7^1 = 0. 

ut 

These last two equations simply state that density 
and viscosity within each fluid remains constant. 

Nondimensionalization gives a Weber and a 
Reynolds number defined by: 

Re = ££E 

O lid 

In addition, the density ratio r = p d l p o and the 
viscosity ratio X = fl d / ji 0 must be specified. 
Here, the subscript d denotes the drop fluid and o 
the ambient fluid. In off-axis collisions, the drops 
approach each other along parallel lines that are 
some distance apart. If this distance is greater than 
the drop diameter, D, the drops never touch and no 
collision takes place. If this distance is zero, we 
have a head-on collision. To describe off-centered 


2 




Figure 1. The computational domain and the initial conditions. The drops 
are initially two and a half diameter apart. 


collision a new nondimensional parameter, usually 
called the impact parameter, is required in addition 
to the Weber and the Reynolds number defined 
earlier. This parameter is usually defined as 



where % is the perpendicular distance between the 
lines that the drops move along before collision. 

The force used to drive the drops together initially 
is taken as 

f x ~ C(p ~ Po){ x ~ x c) 

so the force acts only on the drops. Here C is an 
adjustable constant and x c is midway between the 
drops. This force is turned off before the actual 
collision takes place. Initially, the drops are place 
with their centers two and a half diameter between 
them, and C is varied to give different collision 
velocities. 

To solve the Navier Stokes equations we use a 
fixed, regular, staggered grid and discretize the 
momentum equations using a conservative, second 
order centered difference scheme for the spatial 
variables and an explicit second order time 
integration method. The pressure equation, which is 
non-separable due to the difference in density 
between the drops and the ambient fluid, is solved 
by a Black and Red SOR scheme. Other versions of 
our code use a multigrid iteration. The novelty of 
the scheme is the way the boundary, or the front, 
between the drops and the ambient fluid is tracked. 
The front is represented by separate computational 
points that are moved by interpolating their 
velocity from the grid. These points are connected 
by triangular elements to form a front that is used 
to keep the density and viscosity stratification sharp 
and to calculate surface tension forces. At each time 
step information must be passed between the front 
and the stationary grid. This is done by a method 


similar to the one discussed by Unverdi and 
Tryggvason 12 , that spreads the density jump to the 
grid points next to the front and generates a smooth 
density field that changes from one density to the 
other over two to three grid spaces. While this 
replaces the sharp interface by a slightly smoother 
grid interface, all numerical diffusion is eliminated 
since the grid-field is reconstructed at each step. The 
surface tension forces are computed from the 
geometry of the interface and distributed to the grid 
in the same manner as the density jump. Generally, 
curvature is very sensitive to minor irregularity in 
the interface shape and it is difficult to achieve 
accuracy and robustness at the same time. However, 
by computing the surface tension forces directly by 

F c = crj>t x nds 

we ensure that the net surface tension force is zero, 
or: 

oknds = 0 

Here, n is the outward normal, t a tangent vector 
to the boundary curve for each element and K is 
twice the mean curvature. This is important for 
long time simulations since even small errors can 
lead to a net force that moves the drop in an 
unphysical way. 

As the drops move and deform, it is necessary to 
add and delete points at the front and to modify the 
connectivity of the points, to keep the front 
elements of approximately equal size and as "well 
shaped" as possible. This is described in Unverdi 
and Tryggvason. 12 When the drops are close, we 
rupture the interface, in several of our 
computations, by removing surface elements that 
are nearly parallel and reconnecting the remaining 
ones to form a single surface. Here, this 
restructuring of the interface is done at prescribed 
time if the interfaces are close enough. While this 
rather arbitrary (and we have simply selected the 
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Figure 2. Comparison between a fully three-dimensional simulation (right) and results 
obtained by an axisymmetric code (left). The initial conditions are shown at the top of 
each column and the solution is then shown at three equispaced times for each run. 


restructuring of the interface is done at prescribed 
time if the interfaces are close enough. While this 
rather arbitrary (and we have simply selected the 
time when the drops look close enough) this allows 
some control over the dynamic of the rupture, as 
compared with numerical methods where the front 
is not tracked and the film would always rupture 
once it is thinner than a few grid spaces. For a 
more detailed discussion of this point see Nobari, 
Jan, and Tryggvason 4 . 

The method and the code has been tested in various 
ways, such as by extensive grid refinement studies, 
comparison with other published work and 
analytical solutions. It has also been used to 
investigate a number of other multifluid problems. 
In addition to the computations of head-on 
collisions of drops by Nobari, Jan and 
Tryggvason 4 , Unverdi and Tryggvason 13 simulated 
the collision of fully three dimensional bubbles, 
Ervin 14 investigated the lift of deformable bubbles 
rising in a shear flow (see also Esmaeeli, Ervin, 
and Tryggvason 15 , Jan and Tryggvason 16 examined 
the effect of contaminants on the rise of buoyant 
bubbles and Nobari and Tryggvason 17 followed the 
coalescence of drops of different sizes. Nas and 


Tryggvason 18 presented simulation of thermal 
migration of many two dimensional bubbles. 

Results and Discussions 

For the computations presented here, W£=23, 
Re= 68, r=40, and X =20, but the impact parameter, 
I , is varied. The computational domain is resolved 
by a 32 by 32 by 64 cubic mesh and the drop 
diameter is 0.4 times the shorter dimension. 

While we have done extensive checks of the 
accuracy of our axisymmetric code, the three- 
dimensional code has not been tested as thoroughly. 
We have therefore conducted a few calculations of 
head-on collisions where the results from the three- 
dimensional simulations can be compared with the 
axisymmetric results. Figure 2 shows this 
comparison. The axisymmetric results are to the 
left and the fully three dimensional results to the 
right. The initial conditions are shown at the top of 
each column and the drops are then shown below at 
equispaced times. The force that acts on the drops 
initially is turned off before impact (just before the 
second frame). As the drops collide they become 
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Figure 3. The x-position of the center of mass of 
one drop versus time, as computed by both the 
fully three-dimensional code and an axisymmetric 
one, for two different resolutions. 

flatter, and the ambient fluid between them is 
pushed away, leaving a thin film of fluid between 
the drops. Here, this film is not removed and the 
drops therefore rebound, recovering their spherical 
shape. Obviously, the results are in good 
agreement. Figure 3 shows a more quantitative 
comparison, where we plot the ^-position of the 
center of mass for the drops in figure 2, as well as 
for drops computed on a coarser grid (16 by 16 by 
32). The agreement is reasonably good, although 
the coarse grid results are not in as good agreement 
with each other as the finer grid results are. 

In figure 4, the off-axis collision of two drops, for 
1=0.75, is shown. The pair is shown at several 
equispaced times, beginning with the initial 
position at the top of the figure. Once the drops 
have the desired velocity, around the third frame 
from the top, the force that is applied to drive the 
drops together is turned off. The drops continue to 
move together, and in the fourth frame they have 
collided, deforming as they do so. Since the 
collision parameter is relatively high, the drops 
slide past each other and continue along their 
original path. The bottom four frames show the 
motion of the drops after the collision. During the 
collision the drops become nearly flat where they 
face each other, and as the drops slide past each 
other the fluid layer between the drops becomes 
progressively thinner. If it becomes thin enough it 
should rupture, but here we have not allowed that 
to happen. (As seen in figure 6, rupture of this film 
will change the resulting evolution considerably.) 
In figure 5, the velocity components of the center 
of mass of one of the drops, (a), and the kinetic and 
the surface tension energy, (b), is plotted versus 
time. The solid curve in (a) is the velocity in the 
horizontal direction. It increases as the force 
accelerates the drops together, and then decreases 
slightly due to the drag from the outer fluid after 



Figure 4. Bouncing collision. Here 1=0.75 and the 
drops are not allowed to coalesce. The initial 
conditions are shown at the top and the drops are 
then shown every 0.42 time unit. 
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Figure 5. (a) Center of mass velocity of one drop 
from the computation in figure 4. (b) Kinetic and 
surface tension energy of one drop. 


the force is turned off. When the drops actually 
collide, it is reduced more rapidly, but eventually 
resumes a nearly constant decay rate after the 
collision is over. The velocity component in the 
vertical direction (short dashes) is non zero only 
during the actual collision. The kinetic energy in 
(b) shows similar behavior as the velocity: it 
decreases slowly after the force is turned off, more 
rapidly during collision and then resumes slow 
decay. The surface tension energy rises during the 
collision as the drop deforms, thus contributing to 
the reduction in the kinetic energy. Notice that the 
drop oscillates slightly after the collision as seen in 
the surface tension energy plot 

Although bouncing is observed for real drops, it is 
actually a relatively rare outcome of a collision, 
only seen when the drop deform and trap fluid 
between them and the velocity is sufficiently large 
so the film does not have time to drain before the 
drops rebound. To investigate the behavior of drops 
that coalesce, we have written software to 
automatically remove the front bounding the thin 
film between the drops at a prescribed time and 
allow the drops to coalesce. Figure 6 shows the 
results of two computations where the drops 


coalesce. All parameters are the same as in figure 2, 
except that in the left column 1=0.50, and in the 
right column 1=0.825. The film between the drops 
is ruptured at time 0.46 for both runs. In these 
computations we put t=0.0 when the distance 
between the center of the drops is one diameter. For 
the low impact parameter case, the drops deform 
considerably during the initial impact, as observed 
for head-on collisions, but the impact parameter is 
sufficiently large so the drops still slide past each 
other. As the film is ruptured and the drops coalesce 
the momentum of each drop is sufficiently large so 
the large combined drop continues to elongate. 
Eventually, however, surface tension overcomes the 
stretching and the drop is pulled into a spherical 
shape. Due to the velocity of the drops that 
coalesced, the combined drop rotates. 

While the low impact parameter drops are in many 
way similar to drops undergoing a head-on 
collision, the high impact parameter drops in figure 
6b deform only slightly as they collide. When the 
interface between them is ruptured, they have hearly 
passed each other and after rupture their momentum 
is sufficiently large so they continue along their 
original path and stretch the fluid column 
connecting them until it is near breaking. We have 
not written the software necessary for rupturing the 
filament connecting the drops and therefor must 
stop the computations at this point Notice, that 
the coalesced drop rotates, as the low impact 
parameter one did, although much less. 

In figure 7, the surface tension energy, the kinetic 
energy and the total energy of the drops from figure 
6 are plotted versus time. Initially, the kinetic 
energy is increased by the force that accelerates the 
drops together. Since this force is not constant (it 
increases linearly with distance from the center of 
the computational box) the increase is not quadratic 
as for the computations reported in Nobari, Jan, and 
Tryggvason 11 . After the force has been turned off, 
the drops move a short distance before colliding. 
Since the ambient fluid has a finite viscosity, 
kinetic energy is dissipated due to friction and the 
drops slow down. As the drops come in contact, the 
kinetic energy of the low impact number drops 
decreases rapidly, but the high impact number drops 
are not affected to any significant degree. Similarly, 
the surface tension energy of the low impact 
number drops increases and the drops deform, but 
the surface tension energy of the other drops hardly 
increases at all since the drops remain almost 
spherical. When the film between the drops is 
ruptured, part of the drops surface is removed and 
the surface energy reduced. This reduction is larger 
for the low impact number drops since the area 
removed is larger. Initially, the kinetic energy of 
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Figure 6. Coalescing collisions. The initial conditions are shown at the top of the figure and the pair is then 
shown every 0.42 time units. The film is ruptured at t=0.46 in both cases, but the impact parameter is different 
for the two runs. In the left column 1=0.5, and 1=0.825 in the right one. For the low impact parameter, the drops 
coalesce permanently, but for the higher impact parameter they separate again. 
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1=0.5, (b) 1=0.825. The total energy, the surface 
energy and the kinetic energy of the drops are 
plotted versus time. 

the high impact number drops is nearly unaffected 
(and continues to be dissipated at the same rate as 
before the drops collide), but as the coalesced drop 
starts to stretch and the surface tension energy to 
increase, the kinetic energy drops sharply. As the 
filament between the drops starts to neck down, the 
increase in surface area stops and the kinetic energy 
levels off. For the low impact number drops, the 
rupture takes place near the point of maximum 
deformation and surface energy is initially converted 
into kinetic energy as the drop adjusts to the new 
shape. The momentum of the drop before impact 
is, however, sufficiently large so the drop is 
stretched as the fluid of the original drops continue 
along the paths they were following before 
collision. This leads to an increase in surface 
tension energy and decrease in kinetic energy. When 
the surface tension energy reaches maximum the 
kinetic energy is not zero due to the finite 
rotational motion of the coalesced drop. Eventually, 
the coalesced drop oscillates. 

For modeling of droplet collisions, the major 
question is whether the collision results in a single 


drop or not. Figure 8 shows our computations in 
the I-Re plane. In addition to the computations 
shown in figure 5, we have conducted two other 
calculations , at different impact parameters. The 
runs that lead to a coalesced drops are shown by 
black squares and those leading to grazing collision 
as open squares. We have also plotted the 
experimental results of Jiang et al 3 , for the 
boundary between these two collision modes for 
We=23. Their results do not extend down to the 
Reynolds number simulated here, but since the 
boundary is only weakly dependent on the Reynolds 
number it seems safe to extrapolate their results to 
our Reynolds number. The dashed line shows this 
extrapolation, showing that the numerical results 
are consistent with the experiments. 


Conclusion 

The purpose of this paper is to demonstrate the 
feasibility of accurate numerical predictions of fully 
three-dimensional off-axis collisions of two drops. 
To do so, we have simulated a few cases, both with 
and without rupturing of the interface separating the 
drops. Although the rupture of the film between the 
drops is done in an ad hoc, way, the results are in 
reasonably good agreement with experimental 
observations. For exact predictions of the boundary, 
a more accurate criteria for the rupture time 19 would 
have to be used. These computations, which require 
about ten hours on a CRAY-XMP, are done on a 
relatively coarse mesh and are therefore limited to 
relatively small Reynolds and Weber numbers. 
Nevertheless, they do demonstrate well the 
capability of the method. 
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13. ABSTRACT (Maximum 200 words) 

Three-dimensional simulations of the off-axis collisions of two drops are presented. The full Navier-Stokes equations are 
solved by a Front-Tracking/Finite-Difference method that allows a fully deformable fluid interface and the inclusion of 
surface tension. The drops are accelerated towards each other by a body force that is turned off before the drops collide. 
Depending on whether the interface between the drops is ruptured or not, the drops either bounce or coalesce. For drops 
that coalesce, the impact parameter, which measures how far the drops are off the symmetry line, determines the eventual 
outcome of the collision. For low impact parameters, the drops coalesce permanently, but for higher impact parameters, a 
grazing collision, where the drops coalesce and then stretch apart again is observed. The results are in agreement with 
experimental observations. 
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